///\todo upgrade

//////////////////////////////////////////////////////////////////////////////////
//                                                                              //
//  This file is part of the buola project (https://code.google.com/p/buola/).  //
//                                                                              //
//  Copyright(c) 2007-2012 Xavi Gratal                                          //
//  gratal AT gmail DOT com                                                     //
//                                                                              //
//  Buola is free software: you can redistribute it and/or modify               //
//  it under the terms of the GNU General Public License as published by        //
//  the Free Software Foundation, either version 3 of the License, or           //
//  (at your option) any later version.                                         //
//                                                                              //
//  Buola is distributed in the hope that it will be useful,                    //
//  but WITHOUT ANY WARRANTY; without even the implied warranty of              //
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               //
//  GNU General Public License for more details.                                //
//                                                                              //
//  You should have received a copy of the GNU General Public License           //
//  along with buola.  If not, see <http://www.gnu.org/licenses/>.              //
//                                                                              //
//////////////////////////////////////////////////////////////////////////////////

#ifndef _BUOLA_MAT_SVD_H_
#define _BUOLA_MAT_SVD_H_

#include <buola/mat.h>
#include <buola/mat/detail/eigen_wrapper.h>

namespace buola { namespace mat {

template<typename tDerived>
CMat<detail::OScalar<tDerived>> pseudo_inverse(const CExpr<tDerived> &pExpr)
{
    CMat<detail::OScalar<tDerived>> lResult(pExpr.Derived().Cols(),pExpr.Derived().Rows());
    auto lSVD=detail::eigen_map(detail::eigen_mappable(pExpr.Derived())).jacobiSvd(Eigen::ComputeThinU|Eigen::ComputeThinV);
    typename decltype(lSVD)::SingularValuesType lSV=lSVD.singularValues();
    detail::OScalar<tDerived> lTolerance=1e-6;
    for(int i=0;i<lSV.rows();i++)
    {
        if(lSV(i)<lTolerance)
            lSV(i)=0;
        else
            lSV(i)=detail::OScalar<tDerived>(1.0)/lSV(i);
    }
    detail::eigen_map(lResult)=lSVD.matrixV()*lSV.asDiagonal()*lSVD.matrixU().transpose();
    return lResult;
}

/*namespace mat*/ } /*namespace buola*/ }

#endif
